

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# DISCLAIMER AND GENERAL INFORMATION
#
# File: 1_cbam_main.R
# Purpose: Main analysis
# Date: 10 July 2024
# Data: Survey data (pulled through 0_cbam_prep.R)
#
# Technical disclaimer:
# All analyses in R version 4.4.1 (2024-06-14 ucrt) -- "Race for Your Life"
# R Studio 2024.04.2 Build 764 ("Chocolate Cosmos" Release (e4392fc9, 2024-06-05) for Windows)
# Windows 10 Enterprise, 64-bit
# 12th Gen Intel(R) Core(TM) i7-1255U 1.70 GHz with 16GB RAM
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++




# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (A) Load data and packages ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Run source file to load data
source("0_cbam_prep.R")

# Load packages
library("janitor")
library("ggpubr")
library("tidyverse")
library("modelsummary")
options(modelsummary_format_numeric_latex = "plain")
library("cowplot")
library("gridExtra")

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (B) Descriptive stats ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Basic descriptive information

# Create data set for observations from control group
control.dta <- dta[dta$exp_group_fct=="Control",c("country","cbam_support")]

# Tabulate CBAM support across countries
g <- control.dta %>%
  group_by(country) %>%
  group_map(~tabyl(.x$cbam_support))
names(g) <- unique(dta$country)
g

# Check data is tabulated correctly
table(control.dta$country)
sum(g$germany$n)
sum(g$hungary$n)
sum(g$switzerland$n)
sum(g$uk$n)

# Sample sizes by country (without NAs)
sum(g$germany$n[g$germany$`.x$cbam_support`<6], na.rm=TRUE) # Germany
sum(g$hungary$n[g$hungary$`.x$cbam_support`<6], na.rm=TRUE) # Hungary
sum(g$switzerland$n[g$switzerland$`.x$cbam_support`<6], na.rm=TRUE) # Switzerland
sum(g$uk$n[g$uk$`.x$cbam_support`<6], na.rm=TRUE) # UK

# Country means
control.dta %>%
  group_by(country) %>%
  summarise(mean=round(mean(cbam_support, na.rm=TRUE),2))


# Percentages of approval vs disapproval by country
round(sum(g$germany$percent[g$germany$`.x$cbam_support`<=2], na.rm=TRUE),3)
round(sum(g$germany$percent[(3 < g$germany$`.x$cbam_support`) & (g$germany$`.x$cbam_support`<=5)], na.rm=TRUE),3)

round(sum(g$hungary$percent[g$hungary$`.x$cbam_support`<=2], na.rm=TRUE),3)
round(sum(g$hungary$percent[(3 < g$hungary$`.x$cbam_support`) & (g$hungary$`.x$cbam_support`<=5)], na.rm=TRUE),3)

round(sum(g$switzerland$percent[g$switzerland$`.x$cbam_support`<=2], na.rm=TRUE),3)
round(sum(g$switzerland$percent[(3 < g$switzerland$`.x$cbam_support`) & (g$switzerland$`.x$cbam_support`<=5)], na.rm=TRUE), 3)

round(sum(g$uk$percent[g$uk$`.x$cbam_support`<=2], na.rm=TRUE),3)
round(sum(g$uk$percent[(3 < g$uk$`.x$cbam_support`) & (g$uk$`.x$cbam_support`<=5)], na.rm=TRUE), 3)


# Percentages of indifferent and non-response responses
round(sum(g$germany$percent[c(3,6)]),3)
round(sum(g$hungary$percent[c(3,6)]),3)
round(sum(g$switzerland$percent[c(3,6)]),3)
round(sum(g$uk$percent[c(3,6)]),3)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (C) Figure: CBAM support in control group ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Construct plotting data
df.plot <- as.data.frame(do.call(rbind,g))
colnames(df.plot)[1] <- c("cbam_support")
df.plot$country <- str_sub(rownames(df.plot), end=-3)
df.plot_wona <- df.plot[is.na(df.plot$cbam_support)==FALSE,]

# Color settings: colorblind-friendly palette
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")


# Plot CBAM support for control condition across four countries
q <- ggplot(df.plot_wona) + 
  geom_col(aes(x=factor(cbam_support), y=percent, fill=country, group=country), position="dodge") +
  labs(y="Respondents (percent)", title="CBAM support by country") +
  theme_minimal() +
  theme(legend.position="bottom") +
  scale_x_discrete(name="", breaks=c("1","2","3","4","5"),
                   labels=c("Strongly oppose", "", "Indifferent", "", "Strongly support")) + 
  scale_fill_manual(values=c(cols[1],cols[2],cols[3],cols[4]), name="", labels=c("Germany", "Hungary", "Switzerland", "UK"))
q

# Export plot to PDF (Figure 1)
ggexport(q, filename="./DV_plot.pdf", width=10, height=5)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (D) Experimental results: CBAM support ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Regression models: collapsed treatment conditions
ols_de <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany",]) 
ols_hu <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary",]) 
ols_ch <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland",]) 
ols_uk <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk",]) 

order_full <- c("factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)4"="Treatment: Both")  

models.full <- list(
  "Germany" = ols_de,
  "Hungary" = ols_hu,
  "Switzerland" = ols_ch,
  "UK" = ols_uk)

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- modelplot(rev(models.full), coef_map=order_full, coef_omit="Interc", facet=TRUE) +
  labs(x="Coefficients",
       y="",
       title="Treatment Effects across Countries",
       subtitle="A: CBAM support (1-5 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic() +
  scale_color_manual(values = cols) +
  theme(legend.title=element_blank())
p1


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (E) Experimental results: Employment effect at different levels ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Regression models: collapsed treatment conditions
# Personal impact
ols_de1 <- lm(cbam_person_inv~factor(tr_main), data=dta[dta$country=="germany",]) 
ols_hu1 <- lm(cbam_person_inv~factor(tr_main), data=dta[dta$country=="hungary",]) 
ols_ch1 <- lm(cbam_person_inv~factor(tr_main), data=dta[dta$country=="switzerland",]) 
ols_uk1 <- lm(cbam_person_inv~factor(tr_main), data=dta[dta$country=="uk",]) 


# Regression models: collapsed treatment conditions
# Regional impact
ols_de2 <- lm(cbam_region_inv~factor(tr_main), data=dta[dta$country=="germany",]) 
ols_hu2 <- lm(cbam_region_inv~factor(tr_main), data=dta[dta$country=="hungary",]) 
ols_ch2 <- lm(cbam_region_inv~factor(tr_main), data=dta[dta$country=="switzerland",]) 
ols_uk2 <- lm(cbam_region_inv~factor(tr_main), data=dta[dta$country=="uk",]) 


# Regression models: collapsed treatment conditions
# Country-wide impact
ols_de3 <- lm(cbam_country_inv~factor(tr_main), data=dta[dta$country=="germany",]) 
ols_hu3 <- lm(cbam_country_inv~factor(tr_main), data=dta[dta$country=="hungary",]) 
ols_ch3 <- lm(cbam_country_inv~factor(tr_main), data=dta[dta$country=="switzerland",]) 
ols_uk3 <- lm(cbam_country_inv~factor(tr_main), data=dta[dta$country=="uk",])

order_full <- c("factor(tr_main)2"="Treatment: Jobs")  

models.full <- list(
  "Germany: personal" = ols_de1,
  "Hungary: personal" = ols_hu1,
  "Switzerland: personal" = ols_ch1,
  "UK: personal" = ols_uk1,
  "Germany: regional" = ols_de2,
  "Hungary: regional" = ols_hu2,
  "Switzerland: regional" = ols_ch2,
  "UK: regional" = ols_uk3,
  "Germany: country" = ols_de3,
  "Hungary: country" = ols_hu3,
  "Switzerland: country" = ols_ch3,
  "UK: country" = ols_uk3  
  )  

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p2 <- modelplot(rev(models.full), coef_map=order_full, coef_omit="Interc", facet=TRUE) +
  labs(x="Coefficients",
       y="",
       title="",
       subtitle="B: CBAM effect on person, region, and country (1-5 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic() +
  scale_color_manual(values = cols) +
  theme(legend.title=element_blank())
p2

full <- as_ggplot(grid.arrange(
          arrangeGrob(p1,ncol=2,
            arrangeGrob(p2))
            )
        )

# Export to PDF (Figure 2)
ggexport(full, filename="./Exp_main.pdf", width=11, height=5)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (F) Experimental results: Climate vs trade trade-off ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Basic information about DV
table(dta$tradeoff)

n <- table(dta$tradeoff,dta$country)
table(dta$tradeoff,dta$country) / colSums(n)

table(dta$tradeoff_bin, useNA="ifany")

# Regression models: Collapsed treatment conditions
ols_de <- lm(tradeoff_bin~factor(tr_frame), data=dta[dta$country=="germany",])
ols_hu <- lm(tradeoff_bin~factor(tr_frame), data=dta[dta$country=="hungary",])
ols_ch <- lm(tradeoff_bin~factor(tr_frame), data=dta[dta$country=="switzerland",])
ols_uk <- lm(tradeoff_bin~factor(tr_frame), data=dta[dta$country=="uk",])

# Germany and climate frame
n <- dta %>%
      filter(country=="germany" & tr_frame!=3) %>%
      # group_by(tr_frame) %>%
      # summarise(mean=mean(tradeoff_bin, na.rm=TRUE))
      group_by(exp_group_fct) %>%
      summarise(mean=mean(tradeoff_bin, na.rm=TRUE))
n

# UK and trade frame
n <- dta %>%
  filter(country=="uk" & tr_frame!=2) %>%
  # group_by(tr_frame) %>%
  # summarise(mean=mean(tradeoff_bin, na.rm=TRUE))
  group_by(exp_group_fct) %>%
  summarise(mean=mean(tradeoff_bin, na.rm=TRUE))
n

order_full <- c("factor(tr_frame)2"="Free trade frame",
                "factor(tr_frame)3"="Climate frame")

models.full <- list(
  "Germany" = ols_de,
  "Hungary" = ols_hu,
  "Switzerland" = ols_ch,
  "UK" = ols_uk)

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p3 <- modelplot(rev(models.full), coef_map=rev(order_full), coef_omit="Interc", facet=TRUE) +
  labs(x="Coefficients",
       y="",
       title="Treatment Effects across Countries",
       subtitle="Climate protection versus trade promotion trade-off (0/1 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic() +
  scale_color_manual(values = cols) +
  theme(legend.title=element_blank())
p3

# Export plot to PDF (Figure 3)
ggexport(p3, filename="./Exp_tradeoff.pdf", width=10, height=5)

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (G) Exploratory results: CBAM support by party affiliation ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Correct the vote intention column
dta_vote <- dta %>% 
  mutate(vote_intention=if_else(as.numeric(vote_intention) < 9995,vote_intention, NA)) %>%# Drop don't know etc.
  mutate(vote_intention=if_else(as.numeric(vote_intention) !=7,vote_intention, NA)) %>% #7 is coded as NA
  mutate(vote_intention=if_else(as.numeric(vote_intention) !=80,vote_intention, NA)) %>% # drop free text Hungary vote
  mutate(vote_intention_labelled=as_factor(vote_intention, levels="labels")) %>% 
  mutate(vote_int_corrected=if_else(country=="switzerland", NA, vote_intention)) %>% #delete incorrect Swiss values (rather to referendum intention not vote intention)
  mutate(vote_intention_corrected=case_when(country=="hungary" & vote_intention==10 ~100, #Correct coding of Hungarian data, currently using Swiss recoding
                                            country=="hungary" & vote_intention==20 ~200,
                                            country=="hungary" & vote_intention==30 ~300,
                                            country=="hungary" & vote_intention==40 ~400,
                                            country=="hungary" & vote_intention==50 ~500,
                                            country=="hungary" & vote_intention==60 ~600, 
                                            TRUE ~ vote_int_corrected)) %>% 
  #Recode factor variable with all parties
  mutate(vote_intention=fct_recode(as.factor(vote_intention_corrected),"CDU/CSU" = "1",
                                   "SPD" = "2",
                                   "LINKE" = "3",
                                   "GRÜNE" = "4",
                                   "FDP" = "5",
                                   "AfD" = "6",
                                   "DK-Jobbik-Momentum-MSZP-LMP-Párbeszéd" = "100",
                                   "Normális Élet Pártja" = "200",
                                   "Magyar Kétfarkú Kutya Párt" = "300",
                                   "Megoldás Mozgalom" = "400",
                                   "Mi Hazánk Mozgalom" = "500",
                                   "Fidesz-Magyar Polgári Szövetség-Kereszténydemokrata Néppárt" = "600",
                                   "Conservative Party" = "101",
                                   "Labour Party" = "102",
                                   "Scottish National Party" = "103",
                                   "Liberal Democrats" = "104",
                                   "Reform UK" = "105",
                                   "Plaid Cymru" = "106",
                                   "Green Party" = "107",
                                   "Democratic Unionist Party" = "110",
                                   "Sinn Féin" = "111")) %>%
  select(country, vote_intention, co2_tax_support)


# Add party affiliation data

green <- c("GRÜNE", "Green Party") # GRÜNE, Green Party
main_left <- c("SPD", "Labour Party") #SPD, Labour Party
main_right <- c("CDU/CSU", "FDP", "Conservative Party", "DK-Jobbik-Momentum-MSZP-LMP-Párbeszéd") #CDU/CSU, FDP, conservative Party, DK-Jobbik-Momentum-MSZP-LMP-Párbeszéd
populist_right <- c("AfD", "Fidesz-Magyar Polgári Szövetség-Kereszténydemokrata Néppárt") #AfD, Fidesz-Magyar Polgári Szövetség-Kereszténydemokrata Néppárt

dta_party <- dta_vote %>% 
  mutate(party_type=case_when(vote_intention %in% green~ "Green",
                              vote_intention %in% main_left~ "Main Left",
                              vote_intention %in% main_right~ "Main Right",
                              vote_intention %in% populist_right~ "Populist Right",
                              TRUE ~ NA))

# Create data frame for plotting CBAM support by party affiliation
plot_cbam_data <- dta_party %>% 
  filter(co2_tax_support!=3) %>% # Filter out undecided respondents
  mutate(test_labelled=as_factor(co2_tax_support, levels="labels")) %>% 
  mutate(carbon_law=factor(case_when(co2_tax_support== 4|co2_tax_support== 5 ~ "Support", 
                                     co2_tax_support== 1|co2_tax_support== 2 ~ "Oppose",
                                     TRUE~as.character(co2_tax_support))))  

# Group dataframe and drop NAs
cbam_data <- plot_cbam_data %>% 
  filter(!is.na(party_type)) %>% 
  mutate(country=case_when(country=="germany"~ "Germany",
                           country=="uk"~ "UK",
                           country=="hungary"~ "Hungary")) %>%
  group_by(country, party_type, carbon_law) %>% 
  summarise(count =n()) %>% 
  mutate(perc = round(count/sum(count)*100, digits = 1))

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")


# Create plot of CBAM support by party affiliation (Figure 4)
p1 <- ggplot(data = cbam_data[cbam_data$country!="UK",], aes(x = party_type, y = perc, fill = factor(carbon_law, levels = c("Support", "Oppose")))) +
  geom_bar(stat = "identity", position = "fill", width = 0.6) +
  geom_text(aes(label = paste0(perc, "%; ", "n=", count)), position = position_fill(vjust = 0.5), size = 3.5, color = "black") +
  facet_grid(rows = vars(country), switch = "y") +
  labs(title = "Support for CBAM by Party Affiliation",
       subtitle = "Question: Do you support or oppose the EU Proposal to Introduce an EU Carbon Border Tax?",
       x = "", y = "") +
  scale_fill_manual("", values = c("Support" = cols[2], "Oppose" = cols[1])) +
  theme_classic() +
  theme(legend.position="bottom")
p1 

# Export plot to PDF
ggexport(p1, filename="./CBAM_vote_intention.pdf", width=10, height=10)



# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#                           END HERE
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
